# This file and accompanying functions are renamed with domestic RC added to 
# a CF on BLP data
#
rm(list=ls())
library(tictoc) #tictoc masks data.table::shift() so put it before HeadR (which libraries data.table) 
library(HeadR)
library(doParallel)
library(SQUAREM)
library(fixest)
setwd("~/Dropbox/BLP_REStat/Replication")
library(Rcpp)
# create the Rcpp function we need for computing equilibrium
Rcpp::sourceCpp("Cppfiles/crossDelta.cpp") # gives a warning but don't worry
#parallel processing set up for doParallel, doRNG
cl <- detectCores()-1
registerDoParallel(cores=cl)
getDoParWorkers()
MAXIT <- 500
seedn <- 777 # 666
sequences <- "normal" # alternatives: normal, halton or "sobol"
#
n.hh <- 100000 # careful, this is slow, but works
U0 <- 1  # needed in prob.mat(), but not prob.mat.namedxvars()
n.reps <- 1 # for richsubs (no need to do a lot if a lot of consumers.)
sigmadomesticpct.list <- c(0,3,25,34,50,72,100,125,150,172,200,225,250) # 
num.loops <- length(sigmadomesticpct.list)
#
# This is intended to choose which country is targeted
source(file = "Rfiles/CF_BLPdata_brandiso.R")
target.list <- brands[iso!="USA",unique(iso)] # all non US brands
rm(brands)
#
#
# Here choose which set of parameters to use. 
param.all <- fread("Data/BLP99/pyblp_domestic_param.csv") # estimated using C&G package
param.list <- names(param.all[,par_abbrev:=NULL])
param.list 
#add pub_param (the BLP 99 parameters)
param.chosen <- c("Domestic_pub_param")
#
#
tic()
##
#Loop over sigmadomesticpct list, conducting counterfactuals ====
##
DF <- foreach(i=1:num.loops, .combine = 'rbind') %dopar% {
#
source("Rfiles/CF_BLPdata_domestic_functions.R",local=TRUE)
source("Rfiles/CF_MMX_functions.R",local=TRUE)
source("Rfiles/CF_MMX_MLOG_functions.R",local=TRUE)
#  
sigmadomesticpct <-  sigmadomesticpct.list[i]
cat("Now running for sigma domestic =", sigmadomesticpct ,"\n")
#
# This is used in the counterfactuals : number of characteristics (outside the constant)
if(param.chosen %like% "^Domestic"){
n.x <- 5 #THIS NEEDS TO BE 5 WHEN ADDING DOMESTIC, 4 otherwise 
} else n.x <- 4
#
###
# What is the average own elasticity? ====
###
setwd("~/Dropbox/BLP_REStat/Replication")
source("Rfiles/CF_BLPdata_domestic_getdata.R",local = TRUE) # get original data: original published_param.csv is 1995 and should not be changed (used in other files)
#
year.gph <- 1990
year.monte <- year.gph
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix stuff
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#initialization phase
set.seed(seedn)
data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
BLP.output.m <- BLP.output[m==j]
BLP.meanownelas <- BLP.output.m[, mean(own_elas_m)] # super assignment needed here because of scope for DGP.BLP
BLP.meanownelas
#
###
# Implement CF on BLP ====
###
#
source("Rfiles/CF_BLPdata_domestic_getdata.R",local = TRUE) # get original data
year.gph <- 1990
tar <- 0 # original tariff (to be raised)
tar.increase <- 0.1 # the actual increase in tariff rate
ybar <- exp(y.mn[year ==year.gph]$V2)
# co ownership matrix 
FM <- BLP[year==year.gph,.(m.text=m,f)]
FM[,m := 1:.N]
co.own.mat <- with(FM,build.co(f))
#do a settings file only for the speed 
settings.var <- data.frame(speed =c(0.3,0.3,0.3))
#
# Run the calibration to get xi and mc
year.monte <- year.gph
objective <- "CF"
market <- 1
#
#Run sims for the 2 relevant versions (no logit)
#EHA:
set.seed(seedn)
CF_BLPdata.simulv2 <- CF_BLPdata.simul(vsim=2,cf.meth="EHA")
set.seed(seedn)
CF_BLPdata.simulv3 <- CF_BLPdata.simul(vsim=3,cf.meth="EHA")
#
CF_BLPdata.simulallv <- rbind(CF_BLPdata.simulv2,CF_BLPdata.simulv3) 
CF_BLPdata.table <- CF_BLPdata.simulallv[,list(sigmadomesticpct,v, chg_true_qty,chg_ces_qty)]
#
saveRDS(CF_BLPdata.table,paste0("Tables/CF_BLPdata/SDpar/CF_BLPdata_simul_",param.chosen,sigmadomesticpct,".rds"))
system('echo "Done with one parallel loop" | mail -s "updatefrom Gonzalez" ckhead@gmail.com')
CF_BLPdata.table
}
toc()
#
saveRDS(DF,paste0("Tables/CF_BLPdata/SDpar/CF_BLPdata_simul_",param.chosen,"_all.rds"))
DF <- readRDS(paste0("Tables/CF_BLPdata/SDpar/CF_BLPdata_simul_",param.chosen,"_all.rds"))
